function GenerateSyntheticDWI()

res_pathname = 'C:\Ziang\Dataset\Synthetic\DWI.raw';

width = 64;
height = 64;
depth = 20;

S0 = 1;

GradientOrientations=[
0.1639 0.5115 0.8435;
0.1176 -0.5388 0.8342;
0.5554 0.8278 -0.0797;
-0.4804 0.8719 0.0948;
0.9251 -0.0442 0.3772;
0.7512 -0.0273 -0.6596;
0.1655 -0.0161 0.9861;
0.6129 -0.3427 0.7120;
0.6401 0.2747 0.7175;
-0.3724 -0.3007 0.8780;
-0.3451 0.3167 0.8835;
0.4228 0.7872 0.4489;
0.0441 0.9990 0.0089;
-0.1860 0.8131 0.5515;
0.8702 0.4606 0.1748;
-0.7239 0.5285 0.4434;
-0.2574 -0.8032 0.5372;
0.3515 -0.8292 0.4346;
-0.7680 -0.4705 0.4346;
0.8261 -0.5384 0.1660;
0.9852 -0.0420 -0.1660;
0.1712 0.2576 0.9509;
0.1471 -0.2885 0.9461;
0.3797 -0.4582 0.8037;
0.7994 -0.2011 0.5662;
0.4179 0.4087 0.8114;
0.8135 0.1198 0.5690;
-0.1324 -0.4364 0.8900;
-0.5839 -0.1421 0.7993;
-0.0942 0.4305 0.8977;
-0.5698 0.1788 0.8021;
0.3050 0.6751 0.6717;
0.5085 0.8394 0.1919;
0.3117 0.9494 -0.0368;
-0.2269 0.9724 0.0539;
-0.0115 0.6886 0.7251;
-0.3465 0.8759 0.3360;
0.7411 0.6696 0.0494;
0.9332 0.2164 0.2870;
-0.6260 0.7279 0.2798;
-0.7667 0.2889 0.5734;
-0.0727 -0.6976 0.7127;
-0.4225 -0.8478 0.3207;
0.2438 -0.7111 0.6594;
0.4325 -0.8842 0.1766;
-0.6879 -0.6747 0.2673;
-0.7896 -0.2304 0.5688;
0.6792 -0.7330 0.0370;
0.9103 -0.3028 0.2823;
0.9929 -0.0448 0.1098;
0.9025 -0.0360 -0.4292;
0.4092 -0.1887 0.8928;
0.6587 -0.0357 0.7516;
0.4236 0.1360 0.8957;
-0.1088 -0.1666 0.9801;
-0.3772 0.0084 0.9261;
-0.0945 0.1581 0.9830;
0.2455 0.9391 0.2407;
-0.0747 0.9527 0.2947;
0.1245 0.8414 0.5259;
0.6798 0.6560 0.3279;
0.7940 0.3865 0.4692;
0.5588 0.5583 0.6132;
-0.4784 0.7054 0.5231;
-0.5620 0.4443 0.6976;
-0.2792 0.5940 0.7544;
-0.1585 -0.9475 0.2778;
0.1617 -0.9612 0.2238;
0.0495 -0.8583 0.5108;
-0.5391 -0.6696 0.5108;
-0.5995 -0.4054 0.6901;
-0.3311 -0.5804 0.7440;
0.6192 -0.7190 0.3157;
0.7565 -0.4632 0.4616;
0.5070 -0.6161 0.6028;
0.9755 0.2201 0.0047;
0.9218 0.2253 -0.3158;
0.8613 0.4894 -0.1366;
-0.9524 0.3052 -0.0000;
-0.8986 0.3000 0.3205;
-0.8150 0.5609 0.1459];

S=zeros(size(GradientOrientations,1),width, height, depth);

for z = 1:depth
    for y = 1:height
        for x = 1:width
            for i=1:size(GradientOrientations,1)
                S(i, x, y, z) = S0;
            end
        end
    end
end

for z = 3:depth-2
    for y = 1:height
        for x = 1:width
            % compute directions
            dir1 = [0, 0, 0];
            dir2 = [0, 0, 0];
            
            if x <= 16
                if y >= 16 && y <= 32
                    dir1 = [1, 0, 0];
                end
            else 
                if y < 24 && y > 16
                    dir1 = [1, 0, 0];
                end
                
                dx = x - 16;
                dy = y - 84;
                dis = sqrt(dx * dx + dy * dy);
                if dis <= 60 && dis >= 52
                    dir1 = [-dy, dx, 0];
                end
                
                dx = x;
                dy = y - 32;
                dis = sqrt(dx * dx + dy * dy);
                if dis <= 64 && dis >= 56
                    dir2 = [-dy, dx, 0];
                end
            end
            
            
            n1 = norm(dir1);
            n2 = norm(dir2);
            
            if n1 > 0
                dir1 = dir1 / n1;
            end
            
            if n2 > 0
                dir2 = dir2 / n2;
            end
            
            % simulate DWI signals based on directions
            if n1 == 0 && n2 == 0
                for i=1:size(GradientOrientations,1)
                    S(i, x, y, z) = S0;
                end
            elseif n1 == 0 && n2 ~= 0
                for i=1:size(GradientOrientations,1)
                    S(i, x, y, z) = S0 * SimulateDWMRI(dir2,GradientOrientations(i,:));
                end
            elseif n1 ~= 0 && n2 == 0
                for i=1:size(GradientOrientations,1)
                    S(i, x, y, z) = S0 * SimulateDWMRI(dir1,GradientOrientations(i,:));
                end
            else
                for i=1:size(GradientOrientations,1)
                    S(i, x, y, z) = S0 * (SimulateDWMRI(dir1,GradientOrientations(i,:))+SimulateDWMRI(dir2,GradientOrientations(i,:)))/2;
                end
            end
        end
    end
end


% save dwi
SaveDWI(res_pathname, width, height, depth, S);

end

function SaveDWI(pathname, w, h, d, dwi)
fid = fopen(pathname, 'w');
for z = 1 : d
    for y = 1 : h
        for x = 1 : w
            fwrite(fid, dwi(:, x, y, z), 'float');
        end
    end
end
fclose(fid);
end